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Summary 

Spatially evolving instabilities in a flat-plate 
boundary layer are computed by direct numeri- 
cal simulation (DNS) of the incompressible Navier- 
Stokes equations. In a truncated physical domain, 
a nonstaggered mesh is used for the grid. A 
Chebyshcv-collocation method is used normal to the 
wall, fourth-order finite differences for the pressure 
equation and fourth-order compact differences for the 
momentum equations are used in the streamwise di- 
rection, and a Fourier series is used in the span- 
wise direction. For time stepping, implicit Crank- 
Nicolson and explicit Runge-Kutta schemes are used 
for the time-splitting method. The influence-matrix 
technique is used to solve the pressure equation. At 
the outflow boundary, the buffer-domain technique 
is used to prevent convective wave reflection or up- 
stream propagation of information from the bound- 
ary. Of the techniques available to force transition, 
the present investigation uses approximations from 
linear stability theory (LST) and the newly devel- 
oped parabolized stability equation (PSE) theory for 
inflow forcing. Comparisons are made to (1) validate 
the numerical techniques, (2) determine the effects 
of grid resolution on the downstream evolving flow, 
(3) determine the effects of physical domain trun- 
cation on the disturbance, (4) determine the sensi- 
tivity of the disturbances to changes in the inflow 
forcing, (5) test the outflow boundary condition, and 
(6) test the accuracy of PSE theory. The answers to 
the above objectives would serve as a guide for future 
DNS and PSE studies of more complex problems of 
interest with an a priori knowledge of the preced- 
ing numerical effects. As a note, the present study is 
concerned with unbounded flow transition. Although 
the related problem of bounded flows may be solved 
in a similar manner, the discussion (and references) 
in the present paper are, for the most part, confined 
to unbounded flows. 

Results from the simulations are first compared 
with those of LST with a parallel mean flow used. 
The computed disturbance amplitudes and phases 
are in very good agreement with those of LST (for 
small inflow disturbance amplitudes). Simulations 
are repeated with a nonparallel mean flow. The ex- 
pected increase in growth rate and wavelength shift 
are observed when compared with the parallel mean 
flow case. A comparison is also made between re- 
sults from PSE theory and DNS. A measure of the 
sensitivity of the inflow condition is demonstrated 
with both LST and PSE theory used to approximate 
inflows on “coarse’' and “fine” grids. Very small dif- 
ferences at the inflow are amplified downstream. Al- 
though the DNS numerics are far removed from PSE 


theory, the results agree relatively well. Finally, a 
small-amplitude wave triad is forced at the inflow, 
and simulation results are compared with those of 
LST to verify the accuracy of the three-dimensional 
(3-D) aspect of the code with a known theory. Again, 
very good agreement is found between DNS and LST 
results for the 3-D simulations, and this agreement 
indicates the disturbance amplitudes are sufficiently 
small that nonlinear interactions are negligible. The 
good agreement between DNS and LST results 
verifies that the 3-D aspect of the code is accurate. 

1 Introduction 

For the past century, numerous investigations 
have been conducted in an attempt to predict the 
transition from laminar to turbulent flow in bound- 
ary layers. Most of this effort stems from the in- 
dependent early theoretical accomplishments of Orr 
(refs. 1 and 2) and Sommer fe Id (ref. 3) at the turn 
of the 20th century. Their achievement, based on 
linearized disturbance equations, is a successful ex- 
ample of classical hydrodynamic stability theory and 
is referred to as the Orr-Sornmorfeld equation. It 
was not until some 20 years later that Tollmien 
(ref. 4) was able to solve the Orr-Sommerfeld equa- 
tion, and this solution led to the calculation of a 
critical Reynolds number for the onset of instability. 
On the same subject. Schlichting (ref. 5) computed 
amplification rates of disturbances in the bound- 
ary layer. Part of the first experimental confir- 
mation of the theory was given by Sehubauor and 
Skramstad (refs. 6 and 7), who used a vibrating rib- 
bon to impress a disturbance into the boundary layer 
and hot wires (which were now available) to take 
measurements. With these contributions (and oth- 
ers) spanning some 40 years, theory and experiments 
now agreed on the initial growth of disturbances. 
Today, we have various mathematical and compu- 
tational tools available to solve the Orr-Sommerfeld 
equation. From this equation, much is now under- 
stood concerning boundary-layer disturbances, more' 
commonly referred to as the Tollmien-Schlichting 
(TS) waves. 

Since its origination, stability theory has gained 
wide acceptance and is now a well-established tool 
in the research and engineering community. Further- 
more, it is from stability theory that the first rea- 
sonably comprehensive method for predicting tran- 
sition was derived, the -method by Smith and 
Gambcroni (ref. 8) and Van Ingen 1 . However, the 


1 Van Ingen, J. L.: A Suggested Semi-Empirical Method 

for the Calculation of the Boundary- Layer Transition Region. 
Rep. no. VTH-74, University of Delft (The Netherlands), 1956. 


method is semiempirical and thus requires some fore- 
knowledge of the flow undergoing transition. The 
true physical problem involves disturbances that in- 
teract in a nonlinear manner in later stages of transi- 
tion, and these disturbances are embedded in a grow- 
ing boundary layer. It is apparent that a method, 
which accounts for nonparallel flow and nonlinear 
interactions, is necessary to predict transition. At 
present, such an all-encompassing method of transi- 
tion prediction is beyond our grasp, but progress has 
been made in recent years. 

In the last decade, much excitement has arisen 
because of the strides that have been made in the- 
oretical developments for predicting stages of tran- 
sition beyond the linear growth stage. Stemming in 
part from pioneering attempts at nonlinear theories 
by Benney and Lin (ref. 9) and Craik (ref. 10), Orszag 
and Patera (ref. 11) and Herbert (ref. 12) derived a 
theory, based on Floquet theory, which accounts for 
an experimentally observed three-dimensional (3-D) 
parametric instability. Although the governing equa- 
tions are linearized and a local parallel flow assump- 
tion is made, remarkable agreement is obtained be- 
tween predictions from this new theory arid experi- 
mental results, in particular for the peak-valley split- 
ting mode identified by Klebanoff, Tidstrom, and 
Sargent (ref. 13) and for the peak-valley alignment 
mode observed by Kachanov and Levchenko (ref. 14). 
These are examples of two distinct and different 
routes to transition that are discriminated based on 
the initial disturbance levels. Since its introduc- 
tion in the early 1980’s and subsequent verification 
throughout that decade, the theory for secondary in- 
stabilities is generally accepted and is now widely 
used by the research community as a tool to fur- 
ther understand and predict transition in boundary 
layers. 

More recently, Herbert (ref. 15) and Bertolotti 
(ref. 10) have devised a nonlinear, nonparallel com- 
putational method based on the so-called “parabo- 
lized stability equations” (PSE's). The full bene- 
fits and limitations of this new theory are yet to be 
realized and are explored somewhat in this paper. 
Prior to development of this theory, the only ap- 
proach to solve the nonparallel, nonlinear boundary- 
layer transition problem was by direct numerical 
simulation (DNS), although researchers have had 
some success with asymptotic methods to solve prob- 
lems in the large Reynolds number limit (Smith 
(ref. 17) and Hall and Smith (ref. 18)). To date, 
most studies using DNS have been limited to the 
temporal formulation, in which a spatially periodic 
computational domain travels with the disturbance 
and the temporal evolution of the disturbance is 


computed. This enabled simulations into the later 
stages of transition (Zang and Hussaini (refs. 19 
and 20) and Laurien and Kleiser (ref. 21)), and thus 
provided a data base of qualitative information that, 
however, lacks the physically realistic spatial repre- 
sentation. Spatial DNS provides needed quantitative 
information about transition. But with spatial DNS, 
obstacles exist that have prevented fully carrying out 
such a study. Among these are the realistic spec- 
ification of inflow and outflow conditions and high 
demands on computational resources. Even with to- 
day’s supercomputers, current resources are insuffi- 
cient to fully simulate transition to turbulence in a 
boundary layer in a spatial setting. However, Rai 
and Moin (ref. 22) have demonstrated that the qual- 
itative characteristics of the transition process can be 
captured with todays computers. 

Yet, progress in spatial DNS has been made by, 
among others, Fasel, Rist, and Konzelmann (refs. 23 
to 26) and Spalart (ref. 27) for boundary-layer flow 
and Danabasoglu, Biringen, and Streett (ref. 28) for 
channel flow. To date, results obtained from spatial 
DNS have been compared qualitatively and, with 
some success, quantitatively to results from linear 
stability theory (LST), secondary instability theory, 
and available experiments. For a more complete list 
of accomplishments in transition prediction through 
the use of DNS, refer to the recent review by Kleiser 
and Zang (ref. 29). 

The goal of the present research effort is to intro- 
duce a spatial DNS approach that adequately handles 
outflow problems that can arise, that properly cap- 
tures the flow physics, and that establishes a para- 
metric understanding of DNS. The accomplishment 
of this goal would lead to potential benchmark solu- 
tions for use with future theories. To accomplish this 
goal, confidence in the numerical techniques must be 
established. In this initial study, results from sim- 
ulations of boundary- layer flow over a flat plate are 
compared with those from LST and PSE theory. 

Symbols 

m inflow amplitudes of forced 

disturbances 

B LST matrix 

C\ coefficients for Runge-Kutta march- 

ing, i = 1,2,3 

Ci matrix coefficients for LST, 

i = 0, 1,2, 3, 4 

C z transformed matrix coefficients for 

LST, i = 0,1, 2, 3, 4 
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D 

F 

F r 

m 

G 

H( u) 
ht 
h x 
Inf 

L(u) 

N b 

N b 

N x ,N y ,N z 

P 

P 

P 

Q 

Q 1 

Rx 

i r 


T 

T n (y) 

Tr 

t 

U, V , W 


u 

Uoo 


collocation derivative operator 

right-hand side of pressure equation 

wave frequency 

similarity dependent variable 

right-hand side of eigenvector 
decomposition technique 

momentum equation operator 

time- step size 

streamwise step size 

modified influence matrix 

momentum equation operator 

size of modified influence matrix 

beginning of buffer domain 

number of streamwise, wall-normal, 
and spanwise grid points 

mean- flow pressure component 

disturbance pressure component 

instantaneous pressure 

eigenvector matrix of D 2 -operator 

inverse of Q- matrix 

Reynolds number based on stream- 
wise coordinate 

Reynolds number based on local 
displacement thickness 

Reynolds number based on inflow 
displacement thickness 

attenuation function for buffer- 
domain technique 

parameter for grid stretching 
normal to wall 

period of disturbance, T = 2n/uj 
Chebyshev polynomial of order n 
matrix trace 
time 

mean-flow streamwise, wall-normal, 
and spanwise velocities 

mean-flow velocity vector, 

U = (t/, 1/, W) 

free-stream velocity 


U, V , w 

u 

u 

x r 


x,y,z 


Vinax 

y 

y 

a,l3 


r 

or 

s* 


A 


Aj;, A 


z 


V 

£ 

$ 

p 

n 

U) 

V 


disturbance streamwise, wall- 
normal, and spanwise velocities 

instantaneous velocities, u — U + u 

disturbance velocity vector, 
u = (u, V, u?) 

transpose of fourth-order penta- 
diagonal matrix 

streamwise, wall-normal, and 
spanwise coordinate directions 

physical far-field boundary distance 

similarity variable, y = y \J R x / x 

spectral domain variable, y € [— 1, 1] 

disturbance streamwise and span- 
wise wave numbers 

computational domain 

computational boundary 

local displacement thickness 

boundary-layer thickness at inflow 

displacement thickness at inflow 

eigenvalue matrix of D 2 -operator 

disturbance streamwise and span- 
wise wavelengths 

fluid kinematic viscosity 

temporary variable, £ = dv/dy 

dependent variable for eigenvector 
decomposition technique 

mean-flow stream function 

pressure-like variable 

disturbance normal vorticity 

disturbance frequency 

divergence operator 


Subscripts: 

max maximum 

n gradient normal to boundary 

r tangential component 

oo free-stream conditions 

Superscripts: 

m Runge-Kutta time step 

(m) higher order derivatives 


n 

full-time-step quantities 

o 

with respect to inflow quantity 

T 

matrix transpose 

(V) 

fifth-order derivative 

(vi) 

sixth-order derivative 

* 

Notation: 

with respect to displacement 
thickness 

B-(i 

Bertolotti PSE with six modes 

C-5 

Chang PSE with five modes 

DNS 

direct numerical simulation 

DL 

DNS with LST inflow 

DP 

DNS with PSE inflow 

LST 

linear stability theory 

PSE 

parabolized stability equation 

RK 

Runge-Kutta 

mis 

root- mean-square 

TS 

Tollmien-Schlichting 


A circumflex over a symbol indicates it is a series 
coefficient. 

2 Governing Equations 

The incompressible Navier-Stokes equations are 
solved in the domain shown in figure 1. The stream- 
wise direction is x, the direction normal to the wall is 
//. and the spanwise direction is z. The correspond- 
ing instantaneous velocities are u = (u, v, w) and the 
pressure is p. The momentum equations are given by 

11 , + (u • V)u = -Vp+ Tv 2 u (1) 

l \ 0 

and the continuity equation by 

V • u = 0 (2) 

where subscripts on the dependent variables denote 
partial derivatives with respect to that subscripted 
variable. The equations are nondimensionalized with 
respect to the free-stream velocity U x , the kinematic 
viscosity v , and some length scale at the inflow (say, 
displacement thickness 6*). A Reynolds number can 
then be defined as R* = U^S'/is. The instantaneous 
velocities u and pressure p may be decomposed into 


mean-flow components, U = ({/, K W) and P, and 
fluctuating components, u = (u, ik w) and p : 

u(x,t) = U(x) + u(x,f) and p(x,t.) = P(x) + p(x,£) 

(3) 

where x = (x,y, z). Thus, the flow field is a com- 
posite of mean and unsteady solution components, 
which are determined and computed in the following 
manner. 

2.1 Mean-Flow Component 

The mean boundary-layer flow on a flat plate 
may be described by the boundary-layer equations, 
which are parabolic in the streamwise (x) direction. 
Although a marching algorithm may be used to 
solve the equations for the mean flow ([/, l/), the 
widely used Blasius similarity profile is employed 
for the present study. A detailed description and 
derivation of the mean-flow equations are provided 
in appendix A. 

2.2 Disturbance Component 

The disturbance, or fluctuating, components of 
equations (3) are determined by solving the form 
of the Navier-Stokes equations that results from our 
substituting equations (3) into equations (1) and (2) 
and subtracting out the mean-flow equations. These 
unsteady, nonlinear disturbance equations are 

u, + (u • V)u+ (U • V)u + (u • V)U = -Vp+ —V 2 u 

R* 

(4) 

and the continuity equation is 

V ■ u — 0 (5) 

with boundary conditions 
u = 0 at y = 0 and u — ► 0 as y —+ oc (6) 

Outflow conditions are provided by parabolizing the 
governing equations (4) over a small portion of the 
downstream computational domain. An illustration 
of this is shown in figure 1. This procedure, known as 
the buffer-domain technique, is described in a later 
section. 

Various analytical and numerical techniques are 
now available to introduce a disturbance into the 
boundary layer. For example, Fasel, Rist, and 
Konzelmann (ref. 26) used time-periodic suction and 
blowing, while Krai and Fasel (ref. 30) used heater 
strips. An alternative form of disturbance forcing 
is to introduce a prescribed time-periodic function 
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at the inflow or the free-stream boundary. For the 
present study, the disturbance forcing takes the form 
of eigenfunctions imposed at the inflow boundary. 
Since the emphasis of this study is to verify the nu- 
merical techniques used in the simulations, a con- 
trolled input is required, which may be used by the 
DNS, LST, and PSE codes. 

The inflow condition Uj n is given by the mean flow 
and a disturbance-forcing function, or 

Uj n = U 0 + u° at x = 0 (7) 

where U 0 is the inflow mean component. For the 
present simulations, the disturbances take the form 
of a linear combination of individual functions: 

m=M n=N 

«•= £ £ K. m ■ Ro 

m=—M 7 i=—N 

(8) 

where A° rn represents the 2-D and 3-D disturbance 
amplitudes, which for the flat-plate boundary layer 
are the amplitudes of Tollmien-Schlichting waves. 3 
is a spanwise wave number, and lu is the real distur- 
bance frequency. Time periodicity is assumed, with 
the period T = 2tt/o;. Also, u£ m (y) represents the 
complex eigenfunctions either found from solving the 
Orr-Sommerfeld and Squire equations or obtained 
from a local approximation of the PSE, and the eigen- 
functions are normalized by the maximum stream- 
wise component. Descriptions of LST and PSE the- 
ory and their numerical solution procedures are given 
in appendixes B and C. 

3 Numerical Methods 


fourth-order central finite differences for the pres- 
sure equation arid compact differences for the mo- 
mentum equations are used on the computational 
domain of AT discrete points. At boundary and near- 
boundary nodes, fourth-order differences are used. 
Although nonuniform grids have been implemented 
and tested, the present study involves the use of 
a uniform streamwise mesh. In this section, both 
differencing methods are discussed. 

The objection to, or difficulty with, using higher 
order schemes comes from the required use of addi- 
tional nodes to achieve the higher accuracy. For cen- 
tral differencing, complications may arise at or near 
a boundary, where insufficient nodes are available for 
the differencing. While approximations by forward or 
backward differences may be used, this may reduce 
the overall global order of the scheme or introduce 
numerical instabilities. 

For the standard central, forward, or back- 
ward finite-difference approximations, the function of 
interest is expanded in a Taylor series as 


fn+\ — fn + h- v fn 



h™ 

ml 


Jm) 

Jn 


+ o(h>; i+1 ) 


(9) 


where f n is the function evaluated at node n; h x is 
the step size, uniform for simplicity; and (m) denotes 
the higher order derivatives. Through expansion of 
neighboring nodes in similar series about node 7? and 
combination with equation (9), fourth-order central- 
difference approximations may be found for the first 
and second derivatives f' n and f '' : 


In this section, the following numerical techniques 
required for the spatial simulation are discussed: 
(1) discretization(s) in the streamwise, wall-normal, 
and spanwise directions; (2) time-splitting procedure, 
from which Poisson equations (2-D) or Helmholtz 
equations (3-D) for the pressure are obtained; (3) the 
eigenvector-decomposition method and the influence- 
matrix method, which are employed to solve for 
the pressure; (4) slip-velocity corrections, which are 
introduced because the pressure equation is inviscid 
to ensure that the tangential boundary conditions 
on velocity remain intact; and (5) buffer-domain 
technique, which is used to prevent wave reflections 
at the outflow. 

3.1 Spatial Discretization(s) 

3.1.1 Discretization in the streamwise di- 
rection. In the streamwise direction (^-direction), 


/' = r^-Un-2 - 8/n-l + «/„+! - f„-2 ) + 0(h*) (10) 

+ 16/ n +l — /n+ 2 ) + 0{hl) (11) 

These approximations are used for the interior nodes 
(i.e., those nodes for which the derivative stencil does 
not extend beyond the boundary nodes). For bound- 
ary and near-boundary nodes, fourth-order forward 
and backward differences formed in a similar man- 
ner are used. An explicit form of these forward- 
and backward-difference relations is provided in 
appendix D. 

As the order of the approximation increases, the 
required number of boundary and near-boundary 
relations and the corresponding required number of 
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nodes per derivative stencil increase. Therefore, to 
achieve higher accuracy without involving the use of 
additional neighboring nodes, compact differencing 
is introduced. 

As originally suggested by Kreiss and Oliger 
(ref. 31) and later discussed for fluid dynamics prob- 
lems by Hirsh (ref. 32), first and second derivatives 
for a compact difference may be approximated by 


f = 

J n 


Do 

1 + j'h'Z.D+D- 


fn 


(12) 


and 


whoi 


f" = 

J n 


D+D_ 


, ' + yJ 1 j-D+D- 


fn (13) 


D°fn — TT~(/n + l “ fn-l] 


D+ fn — ^ ifn + 1 “ fn 


(14) 


D-fn — . (fn fn- 1! 
hr 


Through multiplication of equations (12) and (13) 
by the respective denominators, relations for the 
derivatives may be found: 


g/tt-l + 3 fn + g/u+1 - ^“(/n+1 - fn-l) (15) 


and 


1 

12 



1 

12 


f'n+l 


hi 


(f 


n+1 


2/n + fn~\) 

( 16 ) 


These? equations yield tridiagonal systems, provided 
appropriate boundary conditions are applied. The 
approximations are fourth-order accurate and can 
be solved efficiently by LU-decomposition with 
appropriate backward and forward substitutions. 


To make an accuracy comparison between the 
compact-difference (eqs. (15) and (16)) and the 
central-difference (eqs. (10) and (11)) scheme, Taylor 
series expansions are employed. As Hirsh has shown, 
the truncation errors for the compact differences are 


D( f : i ) = -~h i x f {v) and = 

(17) 


while similar error analyses for the central differences 
yield 

E(fn) = -~hU {v) and £(/") = 

(18) 

Although both schemes are fourth-order accurate, 
the compact-difference scheme should lead to more 
accurate approximations as a result of having smaller 
coefficients on the truncation error. 

As yet, no mention has been made about the 
boundary treatment for the compact-difference 
scheme. At the boundaries, Hirsh (ref. 32) used 
a one-sided fourth-order finite difference. Adam 
(ref. 33) suggested additional boundary relations that 
include near- boundary derivatives in the formula- 
tion; yet, the equations retain the tridiagonal nature. 
However, these relations are third-order accurate and 
indicate no additional benefits, compared with direct 
application of high-order one-sided differences. So, 
for the present compact-difference scheme, one-sided 
fourth-order finite- difference boundary conditions are 
used. (See appendix D.) 

Concerning the boundary condition treatment, 
one might choose the second-order boundary condi- 
tions since a numerical instability could be generated, 
in particular, at the inflow with the use of higher 
order approximations. For the present incompress- 
ible spatial DNS, this problem was not encountered 
with fourth-order boundary conditions; but, in at- 
tempts to use fifth-order boundary conditions or a 
sixth-order compact-difference scheme, a numerical 
instability appeared. Recently, this numerical insta- 
bility for the sixth-order methods has been resolved 
through an alternate boundary condition formula- 
tion by Carpenter, Gottlieb, and Abarbanel (ref. 34). 
This new fifth-order boundary treatment has en- 
abled the use of a sixth-order compact-difference 
scheme. Although numerical difficulties surround- 
ing the sixth-order schemes have been resolved, the 
remainder of this study involves the use of fourth- 
order techniques since additional computational ex- 
pense arises solely from the higher order method. In 
a future study, a comparison of the fourth- and sixth- 
order techniques may be undertaken to determine 
if the accuracy gains with the sixth-order method 
outweigh the additional computational expense. 

3.1,2 Discretization in the wall-normal di- 
rection. Normal to the wall (^-direction), Cheby- 
shev series are used to approximate the disturbance 
at Gauss-Lobatto collocation points. A Chebyshev 
series is used since, as Gottlieb and Orszag (ref. 35) 
have shown, it provides good resolution in regions 
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of high gradients (e.g., near boundaries). Properties 
of Chebyshev series with collocation grids are given 
in Gottlieb, Hussaini, and Orszag (ref. 36). A de- 
tailed discussion of spectral method properties and 
their application is provided by Canuto, Hussaini, 
Quarteroni, and Zang (ref. 37). In the present paper, 
only a brief description of the necessary identities is 
provided. 

The Gauss-Lobatto points for Chebyshev series 

are 

Vi = cos(ni/N y ) (i = 0, 1, ■ * • , N y ) (19) 

where N y is the number of domain intervals (or high- 
est degree of Chebyshev polynomials in the series) 
and N y = N y + 1 denotes the number of collocation 
points. Chebyshev polynomials are defined on the 
interval [—1, 1] and are given by 

T n (y) = cos (n cos - * 1 y) (20) 

where n is the order of the Chebyshev polynomial T n . 
A function f(y) may be represented by a Chebyshev 
series at the Gauss-Lobatto points as 

Ny 

f(y) = E a Mv) ( 21 ) 

71=0 

where a n represents the series coefficients. Deriva- 
tives of the function at collocation points may be 
represented by 


dm) 

dy 




( 22 ) 


where repeated indices indicate summation. The 
derivative matrix D is given by 



Cj Vi - y 3 

= 

Vj 

2(1 - y)) 

— 9 


2Ny + l 

Do,o = 

y 

6 


=0,1 ,---,N y ) 

(j = l,2,---,N y - 1) > 

~ D Ny,Ny 

(23) 


where C{ and Cj = 1 for i, j = 1,2,..., Ny — 1 and 
cn = c~kj =2. Higher order derivatives are simply 

l\y 

multiple powers of D, or 

D p = (24) 

where p is the derivative order. 


Since the spectral interpolation function equa- 
tion (20) is defined on [-1, 1] and the physical prob- 
lem of interest has a semi-infinite domain [0, oo] or a 
truncated domain [0,y m ax], a transformation is em- 
ployed. Studies of spectral methods and mapping 
tranformations in unbounded regions have been con- 
ducted by Boyd (ref. 38) and Grosch and Orszag 
(ref. 39). Here an algebraic mapping is used: 


2/max Sp(l + y) 

2 Sp + 2/max(l — y) 


(25a) 


or 


(2 Sp + ymax)y " 2/max <s p 
2/max (^p T y) 


(25b) 


where y G [0,y ma x)> V £ [ — 1 ? l]i 2/max the nor- 
mal distance from the wall to the far- field boundary 
in the truncated domain, and s p controls the grid 
stretching in the direction normal to the wall. As 
a result of the stretching, the normal derivatives in 
equations (23) and (24) are modified as follows: 

D — mD and D 2 = m 2 D 2 + mm'D (26) 

where the metric is defined as m = dy/dy and 
m! — drn/dy . 

3 A. 3 Discretization in the spanwise direc- 
tion . To simulate the evolution of 3-D disturbances, 
the governing equations must be discretized in the 
spanwise direction ( 2 -direction) in addition to the 
stream wise and wall-normal directions. Some ratio- 
nale in the choice of discretization must be used since, 
with this third dimension, the memory requirements 
and epu cost for a simulation can quickly exceed 
current supercomputer capabilities. From boundary- 
layer experiments (refs. 13 and 14) it has been ob- 
served that a distinct periodic structure is evident in 
the spanwise direction. From this observation, span- 
wise periodicity is assumed, and this periodicity al- 
lows for Fourier series representations. With Fourier 
series, spectral accuracy is obtained in the spanwise 
direction and fast Fourier transforms (FFT’s) or sine- 
cosine transforms may be used, either of which allows 
for the fast computing of derivatives. 

In general, a function /(x, y, 2 , t) is represented by 
a Fourier series expansion in the spanwise direction: 


f(x,y,z,t) 


(N z / 2)— l 

E fn(x,y,t)e in ^/^ 

n=-N,/2 


(27) 
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where N z is the number of Fourier modes, X z = 2ir/{3 
is the span wise wavelength, and (i is a specified 
span wise wave number. To compute a derivative, 
equation (27) is first transformed to Fourier space 
by an FFT. The derivative is computed by 

(28) 


omission leads to 


u t _ u m 


where 


= C[ n H m ( u) + 


/nr m 

~^D 2 (u f + u m ) (30) 

rC 0 


H m (u) = L m (u) + 0[ l H m - l (u) 


where ii tl — njj and an inverse transform is used to 
return the computed derivative to physical space. 

Although the full Fourier representation (eq. (27)) 
is correct if the spanwise direction is periodic, more 
cost-efficient derivatives are computed by cosine and 
sine expansions for the special, yet widely used, case 
of symmetry about 2 = 0 (e.g., wave triads and 
some secondary instability calculations). For the 
simulation problem, even functions (i.e., u, v, and p) 
are expanded with cosine series and odd functions 
(i.e.. tr) are expanded with sine series: 

A\ 

{ u * y.zJ) - y, t)c.w(n2-K/\ z z) 

n — U 

(29a) 

and 

v. 

Ih 0 = ^ Wn(.J\ y. t) sin(n27r/A c ^) (29b) 
/!=() 


Equation (27) is used for a spanwise domain of a full 
wavelength (A-), while the use of the symmetry as- 
sumption with equations (29) permits computations 
on half the domain, or a half- wavelength (A^/2). This 
symmetry assumption decreases the computational 
(cpu and memory) requirements by approximately a 
factor of 2. 


and 


L(u) = (U • V)u + (u ■ V)U + (u • V)u - -Lv*u 

R* 0 

Here represents disturbance velocities at the in- 
termediate RK stages, u"' represents velocities at 
previous RK stages (m = 1,2, or 3), u° repre- 
sents velocities at the previous time step, V 2 , — 
d 2 jdx 2 -I- d‘ 2 /dz 2 , and h t is the time-step size. Re- 
call that D is the derivative (eq. (26)). For a full RK 
stage, the momentum equations with the pressure are 

_ r/i C. ni 

= C["H’"(u) + p D 2 (i."'" 1 + u'") - Vp"‘ +I 

(31) 

Subtracting equation (30) from equation (31) leaves 


u 


1 _ ut 


h) n 


c% 1 . 

R: 


-D 2 (u m+1 - u T ) - Vp 7> 


= -V f 


1 


(32) 


where p is an introduced pressure-like quantity. By 
taking the divergence of equation (32) and imposing 
zero divergence of the flow field at each RK stage 
{rn + 1), a pressure- like equation is obtained: 


vV ,l+1 


1 

w 


(V-ut) 


(33) 


3.2 Time-Splitting Procedure 

For the unsteady disturbance equations (4) to (G), 
a time-splitting procedure is used with implicit 
Crank-Nicolson differencing for normal diffusion 
terms and an explicit third-order Runge-Kutta 
method for all remaining terms. The Runge-Kutta 
(RK) scheme, introduced by Williamson (ref. 40), 
was implemented with the Crank-Nicolson method 
for Taylor-Couette flow calculations by Streett and 
Hussaini (ref. 41). This time-splitting procedure con- 
sists of three intermediate RK stages, each stage of 
the following form. 

1 he pressure is omitted from the momentum 
equation (4) for the fractional RK stage, and this 


which is subject to homogeneous Neumann boundary 
conditions. (See ref. 41.) The solution procedure is 
as follows. The intermediate RK velocities are de- 
termined by solving equation (30). The pressure-like 
correction p m+1 is found by solving equation (33). 
Then, the full RK stage velocities u rn+1 are obtained 
from equation (32). Upon solution of the above sys- 
tem three consecutive times, full-time-step velocities 
u n+1 are determined. The RK coefficients and time 
steps are given by 


~£l 

<4 

c-y 


1 

1 

2 

0 _ 

C 2 

c\ 

c i 


9 

4 

1 

2 

-4 

-Cf 

c\ 

Cl J 


32 

L 15 

1 

5 

153 

“IE - J 
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and 


' h) ' 





► = < 

V2 h t 

> (34b) 

t'S 
fl t J 


< 4 hf. j 



where the sum of the three RK time stages equals 
the full time step hf . 

3.3 Eigenvector-Decomposition Method 

To obtain the pressure-like correction p for the 
2-D and 3-D boundary- layer problems, solutions of 
Poisson equations for each time step arc required. 
For 3-D simulations with spanwise periodicity as- 
sumed, the pressure correction is determined in 
transform space, for which the Fourier coefficients 
are solved. In transform space, the Poisson equations 
become Helmholtz equations. In order to solve the 
equations efficiently, a fast elliptic solver is required. 
For this purpose, the tensor-product, or eigenvector 
decomposition, approach is employed. Danabasoglu, 
Biringen, and Streett (ref. 28) used the eigenvector 
decomposition method for the 2-D channel problem. 
The present solver description is for the 3-D simula- 
tion problem. 

The Helmholtz equations in transform space are 
given by 

V'i :v Pn - blih, = F n (35) 

where 3 n represents the spanwise derivative coeffi- 
cients, or wave numbers, of equation (28). The term 
F n represents the transform coefficients of F\ where 
F is the right-hand side of equation (33), and p n gives 
the transform coefficients of the desired solution, p. 

With respect to matrix operations, ( y , x) ordering 
is used below. Discretized in y and x, the Helmholtz 
equations become 

D "pn + pnX^ — /iff pn — F n (36) 

where D 2 is the Chebyshev-collocation operator in 
equation (26) modified to include boundary condi- 
tions, and X 7 is the transpose of the strcamwisc cen- 
tral finite-difference operator, which for the present 
study is fourth-order accurate and leads to a penta- 
diagonal matrix. The matrix D 2 may be decomposed 
into 

D 2 - QAQ" 1 (37) 

where A is a diagonal matrix of eigenvalues and Q 
is the corresponding matrix of eigenvectors of D 2 . A 
new dependent matrix is introduced and defined as 

*n = Q _1 p« (38) 


Substituting equations (37) and (38) into (36), one 
obtains 

A4>„ + <I>„X r - (jfan = G n (39) 

where Gn — Q ^ F n . Equation (39) is used to solve 
for which is then used in equation (38) to solve 
for j pn. Since the coefficient matrix in equation (39) is 
pentadiagonal for fourth-order streamwise discretiza- 
tion, fast back substitutions result. The solution 
is then transformed through inversion to physical 
space. The derivative matrix D, its inverse, and 
matrices Q and Q -1 are mesh-dependent, matrices 
and need to be calculated only once; the same is 
true of the influence matrix, which is described in 
the next section. To reduce the computational cost, 
planes of the computational domain can be sent to 
the solver for vector izat ion. Sending the entire com- 
putational block may be done and leads to a more ef- 
ficient solver, but the resulting memory requirements 
far outweigh the cost savings. 

3.4 Influence-Matrix Method 

Equations (30) to (33) are solved on a nonstag- 
gered grid. An influence-matrix method is employed 
to solve for the pressure. Street t and Hussaini 
(ref. 41) used the method for the Taylor-Couette 
problem, and later Danabasoglu. Biringen, and 
Streett (ref. 28) used the method for the 2-D channel 
flow problem. Instead of solving a Poisson-Neumann 
problem, two Poisson-Dirichlet problems are solved. 

The solution of the following Poisson-Dirichlet 
problem, which is the pressure-like equation, is 
sought: 

V 2 p=F in F (40a) 

p n — 0 on dr (40b) 

where T is the computational domain, dF is the com- 
putational boundary, and p n indicates a derivative 
of the pressure-like quantity normal to the bound- 
ary dF. To accomplish this, a sequence of solutions 
to the following problem is first determined: 

V\/ - 0 in F (41a) 

p l ~ Sj j on dT (41b) 

for each discrete boundary point xj. The Dirac delta 
function 8ij is defined as Sjj = 1 for i — j and 
8j j = 0 for i ± j. Upon computation of the vectors 

of normal gradients p l n at all the boundary points, 
these vectors are then stored in columns to yield a 
matrix that is referred to as the influence matrix, or 


9 



(42) 


subject to the boundary constraint 


Inf = pn li ] 

where N B is the number of boundary points. 

The influence matrix, which is dense, is of order 
N b x N b for 2-D problems and of order N B x N B x N z 
for 3-D problems, and it is dependent on the com- 
putational mesh only. Since the matrix is depen- 
dent on the mesh, it need be calculated only once 
for a given geometry. However, the memory require- 
ments for the influence matrix for a 3-D problem 
can quickly become overbearing and, thus, eliminate 
the possibility of performing simulations into later 
stages of transition. For example, this single ma- 
trix may easily require 70 Mbytes of memory in the 
early stages of transition of a standard 3-D problem 
with N z = 16. However, since the Helmholtz equa- 
tion (35) is solved in Fourier space, where the coef- 
ficients are independent of each other, this memory 
requirement can be alleviated with a small penalty 
of cpu time (fractions of a second). Through se- 
quential reading of the planes N B x N B of the ma- 
trix from disk, the 70-Mbyte requirement dwindles 
to an acceptable 4-Mbyte size that is now indepen- 
dent of the spanwise discretization. In particular, for 
a Cray supercomputer, buffer in(out) commands can 
be used to read (write) data while the program con- 
tinues to execute. Thus, t he overhead cost is virtually 
negligible. 

The composed influence matrix gives the residuals 
of p as a result of the unit boundary condition 
influence, or 

[I.VfIp ~ Residual (43) 

1 he value of one boundary condition is temporarily 
relaxed so that the problem is not overspecified. This 
is done by setting one column of the influence matrix 
to zero, except for the boundary point of interest, 
which is set to unity. The corresponding residual in 
equation (43) is exactly zeroed. 

The Poisson equation with Neumann bound- 
ary conditions is equivalent to the following solu- 
tion of a Poisson problem and a Laplace problem 
(or Helmholtz problems) with Dirichlet boundary 
conditions. First, solve 

V 2 p l = F in T (44a) 

p 1 = 0 on dr (44b) 

Again, compute the gradients normal to the bound- 
ary pf r This gives the influence of the right-hand 
side F on the boundary. Then, solve 

vV' = o in r (45a) 


P — I /VF ' Pn on (45b) 

The final solution that satisfies the original problem 
and the boundary conditions is p = p 1 - p r *. 

Since the gradient, or boundary condition, at one 
discrete boundary point is relaxed in the influence- 
matrix formulation, the desired condition (p n = 0) 
may not hold at that boundary point. In order to 
regain this boundary condition, the pressure prob- 
lem (eqs. (44) and (45)) is resolved, but this time a 
nonzero constant (say 0.01) is added to the right- 
hand side of equations (44). A pressure correc- 
tion p results. The composite solution satisfies the 
boundary conditions at all discrete nodes and con- 
sists, then, of a linear combination of p and p. This 
combination is found by satisfying the following two 
equations: 

a l Pn + a 2pn = 0 on dT{ (46) 

and 

ai + «2 = 1 (47) 

The final pressure correction p m+1 is then given by 

p m+1 = oip+ (1 - a\)p with ai = p n /{p n - p n ) 

(48) 

Upon solution for p m+ 1 , the full RK time-step veloc- 
ities u m+1 are found via equation (32). As a note, 
the corner points are not included in the discretiza- 
tion and are used in the tangential slip-velocity cor- 
rection only. The pressure at the corners is of mi- 
nor significance and interpolations are sufficient to 
compute these pressures. 

3.5 Slip- Velocity Corrections 

The pressure-like correction equation (33) is an 
inviscid calculation and is well posed, provided that 
boundary conditions on the wall-normal component 
of velocity are enforced. At the end of each full RK 
time step, a nonzero tangential velocity component 
may arise at the computational boundary. This is 
referred to as a “slip velocity.” This slip velocity 
may be made small in magnitude, compared wflth the 
RK step size /i" ! , by a proper choice of intermediate 
boundary conditions. The conditions used herein 
were described by Streett and Hussaini (ref. 41), 
based on the w'ork of Fortin, Peyret, and Temam 
(ref. 42). 

The slip velocities on the boundary for equa- 
tion (32) are 

u” t+1 = u ; - h\ n Vp" i+1 
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on <9r (49) 



where r indicates a tangential component on the 
boundary dT. If u* = 0, then u" ,+ 1 = 0(h[ n ). 
Expanding the gradient term of equation (49) into 
a Taylor scries about t = t m , one obtains 

Vp" l+1 = Vp™ + W{Vp”')t + om") 2 } (50) 

Approximate the time derivative term of equa- 
tion (50) by 

(Vp™)f = Vp \~y~ + oim 2 } (5i) 

h't 

and substitute equations (50) and (51) into equa- 
tion (49). The slip velocity is reduced to 0[(/i™) 1 *], 
and the intermediate boundary conditions that result 
are given by 


* , i in 

u T = u HC + h t 


1 + 


V7 , Wl 

vp r - 


K 


rVo' 




(52) 


where upc = 0 f° r a r ighi wall and Upc = n () for an 
inflow condition or for a wall slot condition evaluated 
at the appropriate time in the RK stage. 

3.6 Buffer-Domain Technique 

The buffer-domain technique for effecting a non- 
reflecting outflow boundary treatment was intro- 
duced by Street t and Macaraeg (ref. 43). The 
technique is based on the recognition that, for 
incompressible flow, the ellipticity of the Navier- 
Stokes equations, and thus their potential for up- 
stream feedback, comes from two sources: the vis- 
cous terms and the pressure field. Examination of 
earlier unsuccessful attempts at spatial simulations 
indicated that upstream influence occurs through the 
interaction of these two mechanisms; strong local ve- 
locity perturbations interact with the condition im- 
posed at the outflow boundary to produce a pressure 
pulse that is immediately felt everywhere in the do- 
main, especially at the inflow boundary. Therefore 
both mechanisms for ellipticity have to be treated. 
To deal with the first source of upstream influence, 
the streamwise viscous terms are smoothly reduced 
to zero through multiplication by an appropriate at- 
tenuation function in a “buffer region,” which is ap- 
pended to the end of the computational domain of 
interest. The viscous terms arc unmodified in the do- 
main of interest. To reduce the effect of pressure field 
ellipticity to acceptable levels, the source term of the 
pressure Poisson equation is multiplied by the atten- 
uation function in the buffer domain. This is akin to 


introduction of an artificial compressibility in that re- 
gion and locally decouples the pressure solution from 
the velocity computation in the time-splitting algo- 
rithm. Thus, in effect , the boundary-layer equations, 
which are parabolic and do not require an outflow 
condition, govern the solution at outflow. Finally, 
the advcction terms arc linearized about the imposed 
mean- or base-flow solution in order that the effective 
advection velocity, which governs the direction of dis- 
turbance propagation, is strictly positive at outflow 
even in the presence of large disturbances. 

The attenuation function used in this work is 
similar to that of references 43 and 28: 




u-Nb) m 

(N x -N b )\iJ 


(53) 


where N ) } marks the beginning of the buffer do- 
main and N x marks the outflow- boundary loca- 
tion. For illustration, the buffer-domain region is 
sketched in figure 1. As shown subsequently for the 
current problems, a buffer-domain length of about 
three streamwise wavelengths is adequate to pro- 
vide a smooth enough attenuation function to avoid 
upstream influence. 

The original buffer-domain implementation of ref- 
erence 43 involved a fully spectral discretization, with 
a spectral multidomain being used in the streamwise 
direction as opposed to the high-order finite differ- 
ences used herein and in reference 28. Thus, early 
testing of the buffer-domain method was doin' in an 
even more sensitive setting. Reference 43 shows a 
number of tests of the method in the context of chan- 
nel flow 7 , albeit they were produced with a code that 
had a slight error and produced a small kink in the 
wall vorticity distribution at outflow 7 . Corrected, the 
fully spectral channel-flow 7 simulation code produced 
results that agree with linear stability theory to five 
significant digits in disturbance growth rate. Addi- 
tional unpublished test cases included simulations of 
Poiseuille-Bcnard flow 7 , in which a strongly unstable 
wall temperature condition w 7 as imposed; the temper- 
ature equation w 7 as included in the solution scheme 
w 7 ith the Boussinesq approximation. For this flow 7 , 
the unstable thermal boundary conditions produced 
large recirculation cells, which in some cases had ver- 
tical disturbance velocities three times larger than 
the imposed Poiseuille base-flow 7 centerline velocity. 
These recirculation cells w 7 ere produced by growth of 
the instability (seeded by numerical roundoff error 
of the computer), a process that is knowm to pos- 
sess a global, rather than convective, instability na- 
ture. The lack of upstream influence even in this ex- 
treme test was confirmed by comparison of vorticity 
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distributions across the channel for channel lengths 
that contained between 5 and 10 cells. 

4 Results 

For the present study, the numerical techniques 
employed for the spatial DNS code are systemat- 
ically verified through comparison with the well- 
established LST and the more recently devised PSE 
theory. First, the present solutions from LST are 
compared with previously published results. Next, 
the results from the DNS with a 2-D disturbance 
forcing for both parallel and nonparallel mean flows 
are compared with the LST results. Third, DNS re- 
sults from 2-D disturbance forcing are compared with 
the PSE predictions. The sensitivity of the inflow 
forcing is demonstrated. Finally, a 3-D simulation 
is conducted and discussed. The DNS results from 
wave-triad forcing are compared with LST results for 
small-amplitude disturbances. The authors extend 
their thanks to Gokhan Danabasoglu of the Depart- 
ment of Aerospace Engineering Sciences, University 
of Colorado at Boulder, for the use of his channel sim- 
ulation code. Thanks also go to Fabio Bertolotti at 
the Institute for Computer Applications in Science 
and Engineering, Hampton, Virginia, for supplying 
his PSE results. 

4.1 Solutions of LST 

Although additional documentation of results de- 
rived from LST is arguably unnecessary in this era, 
for completeness and since a comparison of LST re- 
sults with the DNS results is a part of this study, 
a brief independent code verification is performed 
lor solving the Orr-Sommerfeld Squire problem, as 
described in appendix B. 

Results of the Orr-Sommerfeld equation are well 
documented in the literature. Herein, comparisons 
are made with the results of Jordinson (ref. 44), 
who list'd a, finite-difference approach. For a 2-D 
disturbance with Reynolds number R* of 998 and 
frequency u; of 0.1 122, .Jordinson found the stream- 
wise wave number a of 0.308G - 20.0057. If an 
a priori approximation of the eigenvalue is unknown, 
the spectral global method provides an initial esti- 
mate of the eigenvalue. With an initial guess (say, 
0.3086 — A). 0057), the local method is used to re- 
fine the eigenvalue. Convergence results for the lo- 
cal refinement method are shown in table 1. The 
present results are in good agreement with those of 
Jordinson. 

b igure 2 shows the corresponding eigenfunctions 
for the above parameters. Good agreement occurs 


Table 1. Eigenvalues From LST 
[R* = 998; u; = 0.1122; y max = 75; sp = 10] 



o 

32 

0.308C817 - 10.0055527 

36 

0.3086085 - 10.0057926 

40 

0.3086050 - (0.0056964 

44 

0.3085825 - 10.0057164 

48 

0.3085946 - 10.0057069 

52 

0.3084899 - 10.0057088 

56 

0.3085920 - 10.0057083 

60 

0.3085912 - 10.0057084 


in this comparison of eigenvalue and eigenfunction, 
which demonstrates that sound results of LST are 
available for the DNS verification. 

4.2 Comparison of 2-D DNS and LST 

In this section, the accuracy of the numerical 
methods used for the DNS calculations is tested 
for small-amplitude disturbances through compari- 
son with LST results. Initially, a parallel mean flow 
is assumed. Although this is a physically unrealistic 
flow, it adequately mimics the LST assumptions and 
provides a good initial test case. A Reynolds number 
/?* of 900 and wave frequency F r of (u?/. R*)x 10 6 = 86 
at the inflow are chosen somewhat arbitrarily for the 
test case. In an attempt to determine the grid resolu- 
tion requirements, computations are performed on a 
variety of grids from 40/ j; x 41 to 100/ x x 61 (stream- 
wise x wall-normal), where l x refers to the number of 
TS streamwise wavelengths included in the domain 
and 40/ x denotes 40 grid points per wavelength. If, 
for example, l x = 3, then the grid for 40/ x consists of 
120 points in the streamwise direction. The results 
obtained from each grid are in agreement. 

In the physical domain, the streamwise compu- 
tational domain length is varied, depending on the 
number of TS wavelengths of information required. 
Normal to the wall, the domain length is fixed and ex- 
tended from the wall to an upper truncation distance 
where the far-field boundary conditions arc imposed. 
For parallel flow', the far-field boundary is varied from 
y* — 50 to 100 (where y* = y/tf*). A concern with 
the primitive variable formulation lies in the pres- 
sure calculation, w r hich incidentally is avoided w r ith 
the velocity- vorticity approach as a result of not hav- 
ing to solve for the pressure quantity. If the far-field 
boundary is an insufficient distance from the wall, an 
erroneous disturbance arises throughout the compu- 
tational domain. This erroneous disturbance arises 
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as a result of enforcing the far-field boundary condi- 
tions too close to the wall. Similar errors arise with 
different numerical procedures. For example, distur- 
bances in a boundary layer exponentially decay when 
approaching the far field. Using a shooting proce- 
dure, one can integrate the LST equations from the 
wall to the far-field boundary and match the com- 
puted solutions with asymptotically known solutions. 

If this matching is performed an insufficient distance 
from the wall, the computations will not converge 
to the correct solution. Similarly, the present error, 
which can arise from the far-field boundary condi- 
tions imposed an insufficient distance from the wall, 
leads to incorrect results. 

From the computations with a parallel flow, a far- 
field boundary of y * = 50 appears to be the minimum 
distance for an acceptable disturbance error. Normal 
to the wall, grid stretching is used for the boundary- 
layer computations in order to obtain meaningful re- 
sults efficiently. Stretching factors 3p of 6 to 12 are 
chosen to provide a dense distribution of collocation 
points near the wall. (A smaller Sp clusters more 
points near the wall.) The number of time steps per 
TS wave period is varied from 200 to 1000 in order 
to arrive at a rational choice of the time-stop size 
required. Visual agreement of the results is found 
for each of the time-step test cases. ( r l his agreement 
translates to no more than 0.1-percent error.) Since 
the time-splitting procedure is third-order accurate, 
larger time-step sizes may be used (compared with 
those of a second-order Adam-Bashforth method). 
Computations of the present type involve numer- 
ous parameters (e.g., three-directional grid, far-field 
boundary location, streamwise domain length, and 
time-step size). To remove one of these parame- 
ters, a small time-step size is chosen. Hereafter, 320 
time steps per period are used to maintain temporal 
accuracy through the nonlinear simulations. 

As a first example, the streamwise direction con- 
sists of approximately 7 TS wavelengths with 40 grid 
points per wavelength. Further, these seven TS wave- 
lengths are subdivided into a physical domain of 
four wavelengths and a buffer domain of three wave- 
lengths. For the inflow, a 2-D disturbance described 
by equation (8) with amplitude A° q of 0.1 percent is 
forced. The solutions of the Orr-Sommerfeld equa- 
tion are used for the disturbance profile. The in- 
flow forcing is turned on abruptly, and the results of 
the simulation after three and eight periods of forc- 
ing at the inflow boundary are compared with LST 
predictions in figure 3. The computed phase and am- 
plitude for the streamwise u and wall-normal v dis- 
turbance velocity components with downstream dis- 
tance are in agreement with the LST results. After 


eight periods of forcing, the leading wave has exited 
the computational domain without wave reflections. 
This is an indication that the buffer-domain tech- 
nique is functional. From LST, the spatial growth 
rate is a * = -0.004509. Growth rates from the DNS 
are calculated by a simple central- difference approxi- 
mation with the local maximum disturbance stream- 
wise velocity component u mRX ; this simple approxi- 
mation yields the results and errors shown in table 2 
for various grids. Very good agreement is found be- 
tween LST and the present DNS results, compared 
with results from the crude differential method used 
to compute a*. 

Table 2. DNS Growth Rates From Simple 
Central-Difference Approximations 


N x x N u 


Error 0 , percent. 

40/ x x 11 

0.004438 

1 .57 

60 l x x 41 

.004473 

.80 

Tf 

X 

o 

OO 

.004494 

.33 

401 r x 61 

.004440 

1.53 

601j- x 61 

.004473 

.80 

801,- x 61 

.004494 

.33 


“ Error based on comparison with LST growth rate of 
= —0.004509. 

As demonstrated in figure 3, the buffer-domain 
technique has permitted waves to exit the outflow 
boundary without wave reflection. This is accom- 
plished by specifying a buffer domain of three TS 
wavelengths. We determine this length by compar- 
ing the computed results using various buffer re- 
gions with LST. To demonstrate the effects of using 
a buffer domain of insufficient length, the previous 
DNS results of figure 3 arc shown with erroneous 
results in figure 4. The incorrect results occur for 
a buffer-domain length of one TS wavelength. A 
number of buffer-domain parameter variations may 
be found that are adequate to implement the out- 
flow conditions. The length of the buffer domain, 
the number of grid points, and the slope of the at- 
tenuation function are the important elements that 
may be varied. It is likely that having a small slope, 
and a small change in slope of the function relative 
to the grid spacing is of the most importance; how- 
ever, this postulation has not been confirmed by a 
parameter study. With the present attenuation func- 
tion (eq. (53)), the slope is governed by the buffer- 
domain length and becomes smaller with length in- 
crease. Hence, the three- wavelength domain provides 
an adequate outflow region, while the one- wavelength 
domain does not. 


13 




Next, a nonparallel mean flow is used and the 
simulations are repeated. To ensure accuracy, 60 or 
more points per wavelength are used hereafter. For 
the nonparallel mean flow, computations with four 
periods of forcing are conducted. The streamwise u 
and wall-normal v velocity amplitudes of the distur- 
bance, as computed by simulation and with LST, are 
shown in figure 5. The change in the length scale, 
as a result of the growing boundary layer, is evident 
and leads to an increase in growth rate and a shift in 
wavelength. Since R* 0 = 900 and F r = 86 correspond 
to a growing mode near the lower branch of the neu- 
tral curve, increasing growth rates are expected with 
downstream propagation. This is consistent with the 
results in figure 5. Since LST neglects nonparallel 
effects, exact quantitative agreement is not expected 
here. 

In summary, the numerical techniques used for 
spatial DNS were tested by a comparison with LST. 
This comparison was made because LST provides an 
adequate tool to verify the spatial simulation results 
for small-amplitude disturbances. Also, since LST is 
universally accepted and is a well-established theory, 
it lends credence to the DNS results. A paramet- 
ric study was conducted to determine the effects of 
grid refinement and domain size and to determine an 
adequate time-step size. Furthermore, the outflow 
boundary treatment was successfully tested. A com- 
parison with LST is limited in scope because of the 
underlying assumptions of the theory. Better insight 
into the flow physics of transition and a better un- 
derstanding of the DNS numerics could be achieved 
if results from DNS were compared with a more 
complete theory or experiments. 

4.3 Comparison of DNS and PSE Theory 

Recently, a new theory (PSE) has emerged that 
a< counts for boundary-layer growth and nonlinear 
distm banco interactions. In this section, the results 
from spatial DNS are compared with PSE theory 
predictions. First, the effects of inflow disturbance 
variations and grid refinement on the solutions of the 
DNS in the linear and nonlinear regime are discussed. 
Second, DNS results are compared with those of PSE 
theory. Inferences are drawn by comparing DNS 
results to the distorting mean flow results of PSE 
theory. 

As with Bcrtolotti (ref. 16), calculations are made 
with an inflow Reynolds number R* of 688.315, a 
frequency fr r of 86, and a 2-D disturbance forcing 
at the inflow with amplitude q of 0.25 percent 
rms. The inflow corresponds to a streamwise location 
pi ior to branch I of the neutral curve, in a region 
of disturbance decay. With this inflow amplitude, 


the disturbance decays initially until branch I of the 
neutral curve is reached, where the wave then begins 
to grow. The disturbance amplitude grows through 
the region of instability. Farther downstream, after 
passing branch II of the neutral curve and entering 
the region of stability, the wave saturates, or decays. 
The task at hand is to accurately predict the growth 
and decay of this evolving wave. 

4-3.1 DNS parameter variation . How dis- 
turbances are ingested into the boundary layer and 
the effects of this ingestion are topics of the study 
of “receptivity.” (See Reshotko, ref. 45.) For the 
present study, the presence of an ingested distur- 
bance is assumed, and the evolution of that dis- 
turbance with downstream distance is computed; 
however, it is of utmost importance to know and 
understand how small changes in the disturbance 
(amplitude or profile) affect the computed down- 
stream evolution. It is generally accepted that small 
differences in disturbance amplitudes at ingestion 
into the boundary layer lead to varying locations of 
transition. Although these differences in the distur- 
bance may be small at the inflow', they may amplify 
downstream. 

To demonstrate the sensitivity of spatial DNS to 
inflow disturbance variations, two specific DNS com- 
putations are performed with inflows from LST and 
PSE approximations of the Navier-Stokes equations. 
Hereafter, these two simulation cases are referred to 
as DL and DP, respectively. Since PSE theory is an 
integral method, its inflow' condition must be pre- 
scribed by some local approximation, such as w r as 
prescribed and used by Chang, Malik, Erlebacher, 
and Hussaini (ref. 46). The DNS computations are 
performed on five different grid and inflow variations 
and are forced for approximately 28 to 31 TS peri- 
ods. These five cases are shown in table 3. These 
test cases give a variation in inflow, grid resolution, 
and wall-normal domain length. 

Table 3. Direct Numerical Simulation Test Cases 


Case 

A’j- X !\y 

Far field y* nmx 

DL-41 

60 lj- x 41 

r 75 

DP-41 

60 l x x 41 

75 

DL-61 

60 lj x 61 

75 

DP-61 

80 L- x 61 

100 

DP-81 

60/j x 81 

75 


Figure 6 shows the maximum streamwise am- 
plitudes of the Tollmien-Schlichting (fundamental) 
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wave ui, the mean-flow distortion uo, and the first 
harmonic U 2 with downstream distance for LST and 
PSE inflows on the grid 60/ x x41 (DL-41 and DP-41). 
In the regions of small amplitudes (linear), the results 
from both inflows are in good agreement. Farther 
downstream, the wave amplitudes increase to levels 
where nonlinearities are significant. A slight varia- 
tion in the results for the two inflows arises; at the 
saturation of the fundamental wave, the difference 
between the two wave amplitudes is 2.5 percent. The 
streamwise u and wall-normal v velocity components 
for the LST and PSE inflow profiles are given in fig- 
ure 7. The LST profile is almost indistinguishable 
from the PSE profile. One might initially overlook 
the infinitesimal differences, but these differences are 
clearly amplified downstream (fig. 6). This suggests 
that the evolution of a disturbance is very sensitive 
to small changes in that disturbance. Careful consid- 
eration of the disturbance inputs is of utmost impor- 
tance when any computed results are compared with 
those of theory or experiments; otherwise, improper 
conclusions could result. 

For the next comparison, the computation grid is 
refined to 60l x x 61 (DL-61) and 80/^ x 61 (DP-61), 
with corresponding far-field boundaries of 2/ r * iax — 75 
and y ^ = 100. The resulting streamwise ampli- 
tudes of the fundamental wave u\, the mean-flow dis- 
tortion uo, and the first harmonic U 2 are shown with 
downstream distance in figure 8. Similar to the pre- 
vious results (fig. 6), the amplitudes agree in the lin- 
ear regime and a maximum discrepancy appears near 
saturation. Altering the far-field distance and refin- 
ing the streamwise grid leads to insignificant varia- 
tion in visual comparisons of the results. However, 
refining the normal grid from 41 to 61 collocation 
points leads to larger saturation amplitudes. This 
effect indicates that the normal grid may not be ade- 
quate. To obtain a grid- resolved solution, a final test 
case (DP-81) is computed for 81 collocation points 
with a PSE inflow. Results obtained on the vari- 
ous grids with a PSE inflow are shown in figure 9. 
The results indicate that a grid- resolved, or nearly 
grid-resolved, solution has been attained for the in- 
flow disturbance considered. Also, note that a coarse 
grid leads to an underprediction of the saturation 
amplitudes for the fundamental wave, the mean-flow 
distortion, and the first harmonic. 

To obtain the results shown in figure 9, a buffer 
domain of three TS wavelengths and 320 time steps 
per period is used. For the DP-81 case, the compu- 
tations are restarted and permitted to continue until 
the leading wave front has exited the outflow bound- 
ary. This successfully demonstrates that the buffer- 
domain technique is functional for the nonlinear 


calculations. Finally, the computations arc restarted 
using 416 time steps per period to determine if the re- 
sults are time accurate. Visual comparisons of these 
results with those of figure 9 reveal no differences 
with the use of different time-step sizes. So 320 time 
steps per period are sufficient for the present test 
problem. 


4.3.2 Results of DNS and PSE theoi-y . In 
this section, the nonlinear spatial simulation results 
are compared with PSE calculations of Chang et al. 
(ref. 46) and Bertolotti (ref. 16). With the approach 
of Chang et al., a parametric study was conducted. 
It was determined that 100 points normal to the wall, 
a normal distance of 100<5oi and 5 modes of the scries 
given by equation (Cl) lead to sufficiently accurate 
results for the present test problem. Any further re- 
finement of the PSE grid or number of series modes 
leads to no visible change in the results. The stream- 
wise step size was chosen from a comparison with a 
method of multiple-scales solution for a linear dis- 
turbance evolution. For appropriate step sizes good 
visual agreement of the results was found. The re- 
sults of Bertolotti were obtained with six modes of 
the scries from equation (Cl). Hereafter, the Chang 
et al. and Bertolotti PSE cases are referred to as C-5 
and B-6, respectively. Bear in mind that questions 
concerning PSE parameterization have not yet been 
fully answered. This is illustrated in the compari- 
son of the C-5 and B-6 results, shown in figure 10. 
The streamwise amplitudes of the fundamental w r ave 
u\, the mean-flow distortion uo, and the first har- 
monic U 2 are shown with downstream distance. Good 
agreement of the tw r o PSE results is found for the 
fundamental-wave amplitude in both the linear and 
nonlinear regimes, except for a small discrepancy in 
the saturation amplitude, which may be attributed 
to small differences in the inflow disturbance (figs. 6 
and 8). However, a significant unexplained difference 
in the C-5 and B-6 mean-flow distortion quantities 
does appear. Similar to the DNS results (fig. 9), the 
C-5 results capture early evidence of the first har- 
monic between R * = 690 and 900 (barely visible), 
while the B-6 results do not. With these differences 
in the PSE results noted, the converged DNS results 
(DP-81) are compared with the C-5 PSE results be- 
low. The PSE results of Bertolotti have been com- 
pared with DNS results in less detail than that of the 
present paper. (See ref. 15.) Good agreement was 
indicated by the Bertolotti comparison. Unlike the 
present study, wherein the same disturbance ampli- 
tudes and profiles are used, the Bertolotti comparison 
involved a matching of the disturbance amplitudes at 
some downstream location. 
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This C-5 c‘ase was selected arid compared witli the 
DP-81 results since it is the most controlled compar- 
ison: both calculations are forced with the same in- 
flow disturbance. In figure 11, maximum streamwise 
amplitudes of the fundamental wave, the mean-flow 
distortion, and the first harmonic with downstream 
distance from DP-81 are compared with the C-5 re- 
sults. Results for both the fundamental wave and the 
first harmonic are in good quantitative agreement 
throughout the linear and nonlinear regions, while 
some discrepancy occurs with the mean-flow distor- 
tion quantity in the nonlinear region. It may be ad- 
vantageous to view this comparison on a logarithmic 
scale. (See fig. 12.) A comparison of this type sug- 
gests that results for the fundamental wave, the first 
harmonic, and the mean-flow distortion are in better 
agreement than is shown in figure 11. Also, the early 
evidence of the first harmonic between R* = 690 arid 
900 is visually drawn out. It is apparent that finite- 
amplitude differences are suppressed while small dif- 
ferences in near-zero amplitudes are exaggerated as 
a result of the logarithmic scaling. 

To further examine the DP-81 and C-5 results, 
disturbance profiles at two streamwise locations are 
presented in figures 13 to 16. Figures 13 and 14 show 
streamwise components at streamwise locations cor- 
responding to local Reynolds numbers of R* — 1413 
and 1519, respectively. Figures 15 and 16 show wall- 
normal components at the same respective stream- 
wise locations. As shown in figure 11, the first down- 
stream location is midway through the calculation, 
where the mean-flow distortion has a sudden rise, 
and the second is near the fundamental- wave satu- 
ration. The pictured mean-flow distortion, funda- 
mental wave, and first, and second harmonic profiles 
predicted by DP-81 and C-5 are in good qualitative 
as well as quantitative agreement, even in regions of 
high gradients. As before, the exception lies in the 
mean-flow distortion quantity. The disturbance pro- 
file's of the streamwise components reveal that the 
DNS results agree well with the PSE results in re- 
gions of positive influence on the mean flow, while 
in regions of negative influence PSE theory predicts 
stronger distortions than does DNS. From figures 15 
and 16, it is apparent that the wall-normal compo- 
nent of the mean-flow distortion computed by DP-81 
is in agreement with the C-5 results near the wall, 
with the discrepancy increasing with distance from 
the wall. Most likely this discrepancy in the results 
is due to homogeneous Neumann conditions imposed 
in the far- field wall-normal component of the mean- 
flow distortion for PSE theory. Unlike the DNS ap- 
proach, this approach leads to a nonzero wall-normal 
mean-flow component in the far field. As a result, 


the mean flow varies from the Blasius (mean) flow. 
This variation is shown in figure 17 by a compari- 
son of results from the far-field Blasius solution with 
those from the PSE solution. The maximum differ- 
ence in the mean flows occurs near the location of 
wave saturation. 

As a final test, a simulation was repeated to deter- 
mine the effects of Neumann far-field boundary con- 
ditions for the wall-normal component. For compu- 
tational efficiency, the DP-61 case was used since the 
results appear sufficiently converged with 61 colloca- 
tion points. (See fig, 9.) The results are shown in fig- 
ure 18 along with the previous DP-61 and C-5 results. 
Changing the far-field boundary conditions results in 
no apparent variation in the fundamental- wave and 
the first harmonic results; however, the Neumann 
boundary condition affects the mean-flow distortion 
quantity slightly. Larger amplitudes of the mean-flow 
distortion result. The streamwise and wall-normal 
disturbance components for the DP-61 case with the 
Neumann far-field condition arc given with the C-5 
results in figures 19 to 22. These results correspond 
f° R* — 1413 and 1519. A careful comparison of 
the present profiles with the previous results (figs. 13 
to 16), which have homogeneous Dirichlct far-field 
boundary conditions, reveals that better agreement 
between results from DNS and from PSE theory is 
found. Most significantly, the streamwise mean-flow 
distortion profiles with the Neumann boundary con- 
dition are in better agreement.. Only slightly better 
agreement is achieved for the wall-normal mean-flow 
distortion profiles. It appears the? PSE theory with 
the Neumann boundary condition has a strong ef- 
fect on the mean-flow distortion and only a mild to 
negligible effect on the fundamental wave and the 
harmonics. 

It is important to understand the differences in 
the DNS and PSE theory numerical methods to prop- 
erly draw conclusions from the' above comparisons. 
For the PSE theory approach, the disturbance is rep- 
resented by a Fourier series, as described in appen- 
dix C. The expiations are solved in coefficient space, 
where the dependent variables are the Fourier coeffi- 
cients. Boundary conditions are imposed on each co- 
efficient independently. For the zero-order coefficient 
(mean-flow distortion), the boundary- layer equations 
result; thus, the natural far-field boundary condition 
is a homogeneous Neumann condition on the wall- 
normal velocity component. For the fundamental 
wave and the harmonics, the homogeneous Dirichlet 
boundary conditions are the natural physical choice 
in the far field. For the DNS approach, the full- 
disturbance equations are solved and boundary con- 
ditions are imposed on the disturbance. A physically 
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realizable assumption is that the disturbances vanish 
in the far field, or free stream; thus, homogeneous 
Dirichlet boundary conditions on the disturbance are 
a “good” choice. When the DNS results are com- 
pared with the PSE theory results, differences should 
appear as a result of the different boundary condi- 
tions used. It is apparent that these differences are 
small and become most apparent in the mean-flow 
distortion quantities. Some DNS results were ob- 
tained with the Neumann boundary conditions used 
in the far field and are given with PSE theory results 
in figures 18 to 22. Again differences in the results 
are found. Evidently, significant differences remain in 
the boundary condition treatment even though both 
approaches use Neumann conditions. For the DNS 
case, a Neumann disturbance boundary condition is 
enforced, while for PSE theory, a Neumann mean- 
flow distortion component is enforced. This variation 
suggests that, the flow may exit the far-ficld bound- 
ary in the DNS approach with nonzero velocities in 
the fundamental wave and the harmonics, while this 
cannot happen when the PSE theory is used. Ba- 
sically, this difference suggests that the DNS results 
will not be identical to the PSE theory results. 

Another possible explanation for the small dis- 
crepancy in the comparison of DNS with PSE the- 
ory is that as the disturbance grows and reaches fi- 
nite amplitudes, an induced pressure gradient arises, 
which can be calculated by DNS. The PSE theory ap- 
proach assumes negligible streamwise gradients, and 
the boundary-layer equations result for the mean- 
flow component; thus, PSE theory cannot account 
for the existence of the induced streamwise pres- 
sure gradient. This explanation to the discrepancy is 
under consideration and may be explored further. 

4.4 Comparison of 3-D DNS and LST 

To demonstrate the extension to allow for 3-D dis- 
turbances with a Fourier series (eqs. (27) and (29)) 
used in the spanwise direction, a final comparison 
is made between 3-D spatial DNS results and LST 
results for the parallel boundary layer. As in sec- 
tion 4.2, an inflow Reynolds number R* of 900 and a 
wave frequency F r of 86 are used. Computations are 
performed on a mesh 60/ x x41 x5 (streamwise x wall- 
normal x spanwise) involving cosine-sine transforms. 
In the streamwise direction, the computational do- 
main is six TS wavelengths long (three physical and 
three buffer), and each time period is divided into 320 
time steps. At the inflow, a 2-D fundamental wave 
with amplitude A ( { () of 0.01 percent and a pair of 
oblique waves each with amplitude A° _ of 0.01 per- 
cent and spanwise wavelength \ z of 207T are intro- 
duced. The results at spanwise locations of z = 0 and 


Z -\ z j\ after four TS periods of wave-triad forcing 
are given with LST results in figures 23 and 24. Good 
agreement is found for the small amplitudes consid- 
ered. As a result of the good agreement between the 
DNS results and the LST results, one can conclude 
that the disturbance amplitudes are sufficiently small 
that nonlinear interactions are negligible. 

5 Conclusions and Future Directions 

In the present paper, a spatial direct numerical 
simulation (DNS) approach has been introduced for 
two- and three-dimensional (2-D and 3-D) boundary- 
layer transition problems. The numerical techniques 
have been tested by comparison of DNS results with 
results from the linear stability theory (LST) and 
from the newly developed parabolized stability equa- 
tion (PSE) theory. Results of the present study are 
as follows: 

1. Resulting wave amplitudes and phase from 
the DNS are in very good agreement with those 
from LST for 2-D and 3-D small-amplitude 
disturbances. 

2. The influence and effect of small differences at 
the inflow have been demonstrated using LST 
and PSE theory profiles at the inflow. Even very 
small differences in amplitude or profile become 
amplified downstream . 

3. In the comparison of DNS results with those 
of PSE theory, good overall quantitative agree- 
ment is found in the amplitudes and profiles. 
Questions of boundary condition treatment have 
arisen. A difference in the far-field boundary con- 
dition treatment for the PSE theory is identi- 
fied and likely leads to the differing mean-flow 
distortion quantities. For transition prediction, 
where integral quantities are of importance, the 
PSE theory is likely to be a useful tool for the 
engineer. 

Simulation studies of transition on swept wings, 
large-amplitude wave-wave interactions, 3-D suction 
and blowing for generating streamwise vortices, and 
subharmonic forced transition are all underway. Fur- 
ther detailed comparisons of PSE theory with spatial 
DNS for 3-D transitioning flows are also in progress. 
All these ongoing studies arc directed toward quanti- 
fying transitional flows, which previously could only 
be solved for qualitative information. 

NASA Langley Research Center 
Hampton, VA 23665-5225 
June 2, 1992 
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Appendix A 
Mean-Flow Equations 


By substituting the velocities into the boundary-layer 
equations, one arrives at the following equation for 
the similarity profile: 


For a laminar boundary layer, an order of mag- 
nitude analysis yields the importance of each term 
in the Navier-Stokes equations. Prandtl (ref. 47) ob- 
tained the first estimate by neglecting terms of or- 
der 1 /R 2 v and higher. This led to the now- famous 
bound ary- layer equations. For a 2-D, incompressible 
flow, these are 


0U_ 0V_ 

Ox + lh/ 


(Al) 


OU 0U_ _ _0P 1 0 2 U 

Ox Oy Ox R x Oy 2 


(A2) 


subject to boundary conditions 


f"\y) + \f{y)f"(.y) = o (A6) 

with boundary conditions from equations (A3), or 

/(0) = /'(0) = 0 and f”(y — ► c») — ► 1 (A7) 

where a prime indicates d/dy . After equations (A6) 
and (A7) are solved, the resulting mean veloc- 
ity profile components for the boundary layer are 
determined from equations (A5) and are given by 

u = f'(y) and V = 1/2 [yf(y) - f] (A8) 


V (x, 0) — V(x< 0) = 0 and U(x , cx>) — U^x) = 0 

(A3) 


Moreover, the displacement thickness <5* may be 
computed and is given by 


The first significant observation by Prandtl was 
that the normal pressure gradient is negligible and 
the pressure is a known function of x , which is 
assumed to be impressed on the boundary layer 
by the inviscid outer flow. The second item of 
importance is that second derivatives in x have been 
lost in the boundary- layer approximation, the result 
being parabolic equations in x. The equations may 
readily be solved computationally through use of a 
marching algorithm with x as the marching variable. 

One of the most famous and widely used solutions 
to the boundary-layer equations (Al) to (A3) is 
the flat-plate similarity solution obtained by one of 
Prandtrs students, Blasius (ref. 48). For a parallel 
free- stream flow over a flat plate, the free-stream 
velocity U 0 c is constant. A stream function is defined 
in terms of a similarity parameter y by 

= (>'U x x) l/2 f(y) (A4) 


where y = ^R l /~. Corresponding velocities are 
defined by 


U 


dw 

By 


and 



(A5) 


S* = 


py 'QC 1 ™ rOO 

(1 ~U)dy=-— (1 -f')dy 

jo /?y jo 


or 


Jim (y - /) = 1.7207678— - 
rJ R l J 2 


- 1/2 

— = 1.7207678/1* 7 

x 


(A9a) 

(A9b) 


A Reynolds number based on this local length scale 
may be defined as 7?* = U^S* /v. With equa- 

tions (A9), the mean flow (eqs. (A8)) can be con- 
sistently determined on the DNS mesh so that DNS 
result differences due to mean-flow variations are 
prevented. 


Although the parabolic boundary- layer equa- 
tions (Al) and (A2), which describe the mean flow, 
can be solved computationally by a marching al- 
gorithm, it is more convenient to use the similar- 
ity formulation equation (A4) and numerically solve 
the ordinary differential equation (A6). For the 
present problem, a fifth-order, fixed-step Runge- 
Kutta method described by Luther (ref. 49) is em- 
ployed. The solutions are then retained on the 
computational mesh. 
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Appendix B 


Linear Stability Theory 

Governing Equations 

The well-celebrated Orr-Sommerfeld (refs. 1 
and 2) and Squire (ref. 50) equations have led to a 
better understanding of the linear region of transi- 
tion. These equations are derived from a lineariza- 
tion of the Navier-Stokes equations (1) and (2). In 
terms of the normal velocity, the Orr-Sommerfeld 
equation is given by 


— V 2 - t/A - ^ ) V'v 
R * dx dt , 


d 2 U dv 
dy 2 dx 


(Bl) 


with boundary conditions 


d \ > 

u, — — 0 at y = 0 and 
dy 


v, 


dv 

Wy 


0 as y — * oo 

(B2) 


Substituting the normal-mode form equation (B5) 
into equations (Bl) to (B4) yields 

v"" + a{y)v" + b{y)v = 0 (B6) 

fl" + c(y)Q + d(y)v = 0 (B7) 

where 

a(y) = —2 (a 2 + /3 2 ) — iR[aU () (y) — a ;] 
b{y) = (a 2 + 0 2 ) 2 + iR(a 2 + li 2 ){nU 0 {y) - uj\ 

+ iRaU''(y) 

c(y) = -{a 2 + 0 2 ) - iR{aU 0 (y) ~ u>] 
d(y) = -iR0U' o (y) 


with boundary conditions 


For 3-D disturbances, the equation for normal vortic- 
ity D, referred to as the Squire equation, is required 
in addition to the Orr-Somnierfeld equation. The 
Squire equation is 


(* 


d 


d 


— V - U 

dx dt. 


0 


dU_dv 
dy dz 


(B3) 


with boundary conditions 


v, v, ft = 0 at y - 0 and v, v , ft — ► 0 


where a prime indicates d/dy and R — 


as i/ 


R l J 2 . 


(B8) 


By introducing a temporary dependent vari- 
able £ — v f , the derivative boundary conditions 
arc removed. Substituting the derivative matrices 
(eqs. (26)), the Chebyshev series, and the temporary 
variable into equations (BG) to (B8) leads to 


fi = 0 at y — 0 and — > 0 as y —► oo 

(B4) 


- Dijvj = 0 with = 0 (B9) 


Numerical Methods 

For the present study, solutions of the LST equa- 
tions (Bl) to (B4) are required for an inflow condi- 
tion. A global method is outlined to determine the 
discrete spectrum of interest, and a local method is 
presented that may be used to track eigenvalues and 
corresponding eigenvectors efficiently. 

Both 2-D and 3-D disturbances are assumed to be 
travelling waves. A normal-mode form of solution is 
assumed and is given by 

{u.O} = {v,Q}(y)e l ^ ax+fiz ~' j)t ' > + Complex conjugate 

(B5) 

where {v,fr} are the complex eigenvectors, ui is the 
real frequency, 0 is the spanwisc wave number, and 
a = a r + icti is the complex streamwise wave number. 
In stability theory, a, gives a measure of the distur- 
bance growth, or decay. The streamwise wavelength 
is defined by A^ = 27 t /a r - For 3-D instabilities, the 
spanwise wavelength is defined by A- = 2n/0. 


and 

( D; 3 , + Dj j (ij ) + bj Vj = 0 with t’o = viy = £<> = 0 

(BIO) 

( l)f . + c^flj + djVj = 0 with fi () = Qjv = 0 
^ ' ' (B 1 1 ) 

where 

( 1 / = — [2(n 2 + d 2 ) + iR{nUj - u>)]I 

bj = [(a 2 + d 2 ) 2 + iaRUj + iR{n 2 + (0)(nUj - u>)]I 

Cj = — [2(o 2 + a 2 ) + iR(nUj - w)]I 

d 3 = -(ianu’j) i 

and I is the identity matrix. 

The spatial stability of the boundary layer is of 
interest. The Reynolds number R, frequency u, and 
spanwise wave number 0 are specified, and the com- 
plex streamwise wave number a is the eigenvalue. 
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Klim mating the dependent variable £ and combin- 
ing equations (B9) to (Bll) results in the following 
matrix eigenvalue problem: 

B (a){v n M n } T = 0 (B12a) 

where v n and are the coefficients of the series, or 
discrete functional values of the normal velocity and 
vorticity at the Gauss-Lobatto points, and 

B(rv) — C jo 4 -f- + Cjq C 0 

(B12b) 

I he matrix coefficients C* are complex square matri- 
ces of order N for 2-D instabilities and of order 2N 
for 3-D instabilities. Matrices C4, C3, C2, and Ci are 
singular. The eigenvalue problem is nonlinear and of 
order four in the eigenvalue a. Various methods are 
available to solve such problems. Four approaches 
are given in some detail by Joslin (ref. 51). Herein, 
a global and a local method are used to generate an 
inflow disturbance forcing for the simulations and are 
described next. 

Global Method 

A global method gives the discrete spectrum of 
eigenvalues without a priori knowledge of the value. 
A method referred to as the linear companion ma- 
trix method was given by Gohberg, Lancaster, and 
Rodman (ref. 52). The method has been applied to 
the 2-D Orr-Sommerfeld problem with a flat-plate 
boundary layer by Bridges and Morris (ref. 53), 
among others. 

The linear companion matrix method is a lin- 
earization of the nonlinear problem. An algebraic 
eigenvalue transformation A = l/(a — ,$), where 
s ~ a;/ 0.35, is somewhat arbitrarily used to remove 
the singularities in the coefficient matrices. The 
linearization yields 


where I is the identity matrix of order 2 N and I 
is the identity matrix of order SN. The matrices 
Ci to C 4 are nonsingular as a result of the applied 
eigenvalue shift. The eigenvalues and corresponding 
eigenvectors are found from equation (B13) by using 
the QR algorithm. 


Local Method 

The second solver is a more efficient local eigen- 
value refiner referred to as the Lancaster refinement 
method (ref. 54). The method requires a sufficiently 
accurate initial estimate of the eigenvalue, which can 
be obtained from the above global method. The 
iterative formula is given by 

a ?+1 - a* - 2/(a i )/[/V) - / (1) (a*)] (B14a) 

wdiere 

f(a') = TYIB- VJB^V)] (B14b) 

/ (1 V) = TY{B-V)B( 2 V) 

- [B- VJB^V)] 2 } (B14c) 

and Tr is the matrix trace, B -i is the inverse of B 
(from eqs. (B12)), and B^ denotes the jth derivative 
of B with respect to a. Upon convergence on the 
eigenvalue, the eigenvector may efficiently be found 
by the inverse iteration formula 

B(a){n* +1 , = a {vK U k } T (B15) 


■- c, 'C, 
I 
0 
0 


-Cj-'Ca 

0 

I 

0 


-c- l c } 

0 

0 

I 


- c.r 1 Co “ 
0 
0 
0 


- AI = 0 
(B13) 


where a is a normalizing factor. The procedure con- 
verges in two or three iterations for an initial guess of 
{t>°, - [1,1 ,..., l] r . Equation (B13) or (B15), 

with the continuity equation and the definition of 
normal vorticity, leads to the eigenfunctions {a, v, w} 
required for the inflow condition (eqs. (7) and (8)). 
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Appendix C 
PSE Theory 

For the present paper, it is important to test and 
to verify the accuracy of the numerical techniques for 
the spatial DNS procedure and to make a detailed 
comparison with results of PSE theory. For this 
reason, brief highlights of some of the important 
theoretical and computational aspects of the PSE 
theory are given below. 

As discussed in section 2, the evolution of dis- 
turbances is governed by the unsteady partial differ- 
ential equations (4) to (6). Instead of solving these 
equations directly as in the DNS approach, the PSE 
theory seeks approximate solutions of the parabo- 
lized version of equations (4) to (6). The approxima- 
tion needed to parabolize the governing equations, 
as first suggested by Herbert (ref. 15) and Bertolotti 
(ref. 16), includes the following two assumptions: 
(1) the dependence of the convected disturbance on 
downstream events is negligible and (2) no rapid 
streamwise variation (i.e., S 2 /<9x 2 « 1) occurs in 
the wavelength, the growth rate, the mean velocity 
profile, and all disturbance profiles. 

For nonlinear disturbances present in the flow 
field, periodicity in both the time domain and 
the spanwise domain is assumed and the total 
disturbance in the following Fourier series expansion 
is 

N z N t 

u(x, y, z,t) = ^ ^ Um,n(ah y)e mi ^ z 

■m=—N z n=-Nt 

(Cl) 

where N z and Nt are total numbers of modes kept in 
the truncated series and uj and (3 are the correspond- 
ing frequency and spanwise wave number. Equa- 
tions (Cl) are for velocity components; a similar ex- 
pansion can be written for the pressure p . Through 
substitution of equations (Cl) into the governing 
equations (4) to (6), a set of elliptic equations for the 
transformed variable u m . n Pm,n is obtained. Because 
of the wave nature of these transformed variables, 
they are decomposed into a fast-oscillatory wave part 
and a slow- varying shape function part as 

Pm,n} = {u m ,n>Pm,n} e m,n (C2) 

The governing equations now reduce to a set of 
partial differential equations for shape functions 


{u mn .p m n }. In equation (C2), the fast-scale vari- 
ation along the streamwise direction x is now rep- 
resented by the streamwise wave number a m , n , and 
therefore the second-order variation of shape function 
in x is negligible (based on assumption 2 above). This 
observation leads to the parabolized stability equa- 
tions for the shape functions, which are obtained by 
neglecting all second derivatives in the streamwise 
direction and the terms associated with upstream in- 
fluence. In other words, through proper choice of 
a miU , the evolution of disturbances can then be de- 
scribed by the parabolized equations for the shape 
functions. 

Based on decomposition equations (Cl) and (C2), 
the linear PSE can be derived for any disturbance 
with given frequency and spanwise wave number . For 
nonlinear problems, the following nonlinear terms 
must be added to the governing equations: 

F(x ,y,z,t) = (u ■ V)u (C3) 

Since in the PSE approach the governing equations 
are solved in the wave number space, equation (C3) 
is expanded to a truncated Fourier series in the wave 
number space. The Fourier coefficients then provide 
a nonlinear forcing to each of the linearized shape 
function equations. These inhomogeneous equations 
for the shape functions are solved by a marching pro- 
cedure along the streamwise direction for all Fourier 
modes. 

Numerically, a second-order backward differenc- 
ing is employed to integrate the equations in the 
streamwise direction. High-order finite-difference 
schemes (fourth-order) are employed to discretize the 
normal derivatives. The form of boundary condi- 
tions required for the PSE approach is of particu- 
lar interest. Similar to the DNS approach, no-slip 
conditions arc applied at the wall. The fundamental 
wave and harmonics vanish in the far field. To ac- 
count for the change of displacement thickness in the 
perturbed boundary- layer flow, the far-field normal 
velocity gradients vanish for the mean-flow distortion 
equations. 

Chang et al. (ref. 46) have extended the PSE 
numerical approach to the study of compressible 
boundary layers. With the PSE approach of Chang 
et al. for M oc = 0 and the incompressible results of 
Bertolotti (ref. 16), the present comparison with the 
DNS results is made. 
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Appendix D 

Finite Difference Relations 

The fourth-order finite-difference derivatives at 
the boundary (0, N x ) and near-boundary ( 1 . N x — 1) 
nodes are listed below. The first derivatives are 

/(> = ( _ 25/o + 48/ 1 - 36/2 + I6/3 - 3/4) 

f i ~ _ 10/i + I8/2 - 6/3 -I- fa) 

/.V — j 2 J 1 ( 25 /,V — 48/,v_i + 56 /y 2 

- 16/jv-a + 3/.v_4) 

/V~ 1 = i2h r + ^ ^/;V- 1 _ f 8/,\ r _2 

+ 6/,v_ 3 - /v-l) (Dl) 


The second derivatives are 

■A) ~ y^2 ( 35 /o _ 104 + 114 ^2 - 56/3 + H/ 3 ) 

= Y^2 ( n /« “ 2 °/l + 5 4 /2 + 4/3 - fa) 

fN — J 2 i^2^(35/v — 104/y_j -)- 114/y_2 
- 56/ jV _3 + H/.V— 4 ) 

//V-l ~ ^2 /j2 (^4/v — 20/,v_i + 64/jy_2 

+ 4/jV — 3 - /jV- 4 ) (D2) 
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Figure 1. Computational domain for transitional boundary-layer problem 








Figure 3. Amplitudes of streamwise and wall-normal velocity components of 2-D disturbance with downstream 
distance for parallel boundary layer for inflow. R* 0 = 900; F r — 86; A\ 0 — 0.1 percent. 



Figure 4. Amplitudes of streamwise and wall-normal velocity components of 2-D disturbance with downstream 
distance for parallel boundary layer for inflow with adequate and inadequate buffer regions. R* 0 = 900; 
F r = 86; = 0.1 percent. 
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Figure 5. Amplitudes of streamwise and wall-normal velocity components of 2-D disturbance with downstream 
distance for nonparallel boundary layer for inflow. R* = 900; F r = 86; A ( [ 0 = 0.1 percent. 
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Figure 6. Maximum streamwise amplitudes of fundamental wave mean-flow distortion uq, and first 
harmonic u 2 with Reynolds number (or downstream distance) for cases DP-41 and DL-41 with inflow. 
R* — 688.315; F r = 86; A° 0 = 0.25 percent. 


28 





percent rms 
.125 




Figure 7. Streamwise and wall-normal velocity components of 2-D disturbance as predicted by LST and PSE 
for inflow. Rl = 688.315; F r = 86; A? n = 0.25 percent. 
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Figure 8. Maximum streamwise amplitudes of fundamental wave uj , mean-flow distortion uq, and first 
harmonic u 2 with Reynolds number (or downstream distance) for cases DP-61 and DL-61 with inflow. 
R*> = 688.315; F r = 86; A " q = 0.25 percent. 
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Figure 9. Maximum streamwise amplitudes of fundamental wave u,, mean-flow distortion u 0 , and first 
harmonic u 2 with Reynolds number (or downstream distance) for cases DP-81, DP-61, and DP-41 with 
inflow. R* = 688.315; F r = 86; A? 0 = 0.25 percent. 






R* 



R* 


Figure 10. Maximum streamwise amplitudes of fundamental wave iti, mean-flow distortion u 0) and first 
harmonic u 2 with Reynolds number (or downstream distance) for PSE cases C-5 and B-6 with inflow. 
R*! = 688.315; F r = 86; A" () = 0.25 percent. 
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Figure 11. Maximum streamwise amplitudes of fundamental wave u i . mean-flow distortion U(>, and first 
harmonic U 2 with Reynolds number (or downstream distance) for cases DP-81 and C-5 with inflow. 
R* = 688.315; F r = 86; >1? n = 0.25 percent. 







Figure 12. Logarithmic maximum streamwise amplitudes of fundamental wave tti, mean-flow distortion uq, and 
first harmonic u 2 with Reynolds number (or downstream distance) for cases DP-81 and C-5 with inflow. 
R* = 688.315; F r = 86; A% ft = 0.25 percent. 
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Figure 13. Streamwise disturbance profiles of mean-flow 
second harmonics uo and U 3 with normal distance fr< 
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Figure 14. Streamwise disturbance profiles of mean-flow distortion no, fundamental wave u\, and first and 
second harmonics u 2 and a 3 with normal distance from wall as predicted by cases DP-81 and C-5 for 
R* = 1519. 
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Figure 15. Wall-normal disturbance profiles of mean-flow distortion u 0 , fundamental wave v\, and first and 
second harmonics Vo and i '3 with normal distance from wall as predicted by cases DP-81 and C-5 for 
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Figure 1G. Wall-normal disturbance profiles of mean-flow distortion , H) . fundamental wave m, and first and 
second Jiarmomcs v 2 and V3 with normal distance from wall as predicted by cases DP -81 and C-5 for 
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Figure 18. Maximum streamwisc amplitude of fundamental wave u \ . mean-flow distortion wq, and first 
harmonic a 2 with Reynolds number (or downstream distance) for cases DP-61 and C-5 with inflow. 
R* = 688.315; F r = 86; .4? () - 0.25 percent. 
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Figure 19. Streamwise disturbance profiles of mean-flow distortion no, fundamental wave ui, and first and 
second harmonics U 2 and U 3 with normal distance from wall as predicted by cases DP-61 ( D ?; l30 = 0) and 
C-5 for R* = 1413. 







percent rms percent rms 


Figure 20. Streamwise disturbance profiles of mean- flow distortion uq, fundamental wave u \ , and first and 
second harmonics i/o and 1 x 3 with normal distance from wall as predicted by cases DP-61 (Dv 0 0 = 0) and 
C-5 for R* = 1519. 
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Figure 21. Wall-normal disturbance profiles of mean-flow distortion i\). fundamental wave v\, and first and 
second harmonics V 2 and (,’3 with normal distance from wall as predicted by cases DP-61 (Dv 0 c = 0) and 
C-5 for R* = 1413. 
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Figure 22. Wall-normal disturbance profiles of n 
second harmonics V 2 and r : ;j with normal dist 
C-5 for R* = 1519. 










Figure 23. Streamwise and wall-normal velocity components of wave-triad disturbance with downstream 
distance at 2 = 0 for 3-D parallel boundary layer for inflow. R* = 900; F r = 86; A r { 0 = 0.01 percent; 
A° ±1 = 0.01 percent; X z = 207r. 
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’lguit' -I. Stream wise, wall-normal, and spanwise velocity components of wave-triad disturbance with 
downstream distance at z = X z /4 for 3-D parallel boundary layer for inflow. R* = 900- F — 86" 
~ 0 01 percent: A‘{ ±1 = 0.01 percent; A s = 20? r. 
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